function [ Fgrid ] = Matrix_Bivariate_4pointRefinement( grid , omega )
%-----------------------------------------------------------------
% Input:
% grid(m,n,N,N) - a grid over Z^2 of N X N matrices
% Output:
% Fgrid(2*m-5,2*n-5,:,:)
%-----------------------------------------------------------------
% ABSTRACT
% Tensor product. the middle point is insert in the order y->x of
% calculation
%-----------------------------------------------------------------
% NIr Sharon, 21-05-12
%-----------------------------------------------------------------

%-----------------------------------------------------------
% Initializing

m = size(grid,1);
n = size(grid,2);
N = size(grid,3);

Fgrid = zeros(2*m-5,2*n-5,N,N);

%CheckArray = zeros(2*m-5,2*n-5);

Firstline = zeros(1,2*n-5,N,N);
Lastline  = zeros(1,2*n-5,N,N);
Firstrow = zeros(1,2*n-5,N,N);
Lastrow  = zeros(1,2*n-5,N,N);

%-----------------------------------------------------------
% First main loops: For vertices (new and old) on the grid
for x=2:(m-2)
    for y=2:(n-2)
        Fgrid(2*x-3,2*y-3,:,:) = squeeze(grid(x,y,:,:));
       % CheckArray(2*x-3,2*y-3) = CheckArray(2*x-3,2*y-3) + 1;
        Fgrid(2*x-2,2*y-3,:,:) = Matrix_FourPointAvarage(squeeze(grid(x-1,y,:,:)) , squeeze(grid(x,y,:,:)) , squeeze(grid(x+1,y,:,:)) , squeeze(grid(x+2,y,:,:)),omega);
       % CheckArray(2*x-2,2*y-3) = CheckArray(2*x-2,2*y-3) + 1;
        Fgrid(2*x-3,2*y-2,:,:) = Matrix_FourPointAvarage(squeeze(grid(x,y-1,:,:)) , squeeze(grid(x,y,:,:)) , squeeze(grid(x,y+1,:,:)) , squeeze(grid(x,y+2,:,:)),omega);
       % CheckArray(2*x-3,2*y-2) = CheckArray(2*x-3,2*y-2) + 1;    
    end
end
% The last vertix:
Fgrid(2*m-5,2*n-5,:,:) = squeeze(grid(m-1,n-1,:,:));
%CheckArray(2*m-5,2*n-5) = CheckArray(2*m-5,2*n-5) + 1;
%-----------------------------------------------------------


%-----------------------------------------------------------
% Bounderies handeling: initializing the first and last lines
for y=2:(n-2)
    Firstline(1,2*y-3,:,:) = squeeze(grid(1,y,:,:));    
    Firstline(1,2*y-2,:,:) = Matrix_FourPointAvarage(squeeze(grid(1,y-1,:,:)) , squeeze(grid(1,y,:,:)) , squeeze(grid(1,y+1,:,:)) , squeeze(grid(1,y+2,:,:)),omega);
 
    Lastline(1,2*y-3,:,:) = squeeze(grid(m,y,:,:));
    Lastline(1,2*y-2,:,:) = Matrix_FourPointAvarage(squeeze(grid(m,y-1,:,:)) , squeeze(grid(m,y,:,:)) , squeeze(grid(m,y+1,:,:)) , squeeze(grid(m,y+2,:,:)),omega);
end    
%-----------------------------------------------------------
% Bounderies handeling: initializing the first and last rows
for y=2:(m-2)
    Firstrow(1,2*y-3,:,:) = squeeze(grid(x,1,:,:));    
    Firstrow(1,2*y-2,:,:) = Matrix_FourPointAvarage(squeeze(grid(x-1,1,:,:)) , squeeze(grid(x,1,:,:)) , squeeze(grid(x+1,1,:,:)) , squeeze(grid(x+2,1,:,:)),omega);
 
    Lastrow(1,2*y-3,:,:) = squeeze(grid(m,y,:,:));
    Lastrow(1,2*y-2,:,:) = Matrix_FourPointAvarage(squeeze(grid(x-1,n,:,:)) , squeeze(grid(x,n,:,:)) , squeeze(grid(x+1,n,:,:)) , squeeze(grid(x+2,n,:,:)),omega);
end    
%-----------------------------------------------------------
% Bounderies handeling: vertices depends on one direction alone
% Here - 'y' fixed
y = 2*n-5;
for x=2:(m-2)
        Fgrid(2*x-3,y,:,:) = squeeze(grid(x,n-1,:,:));
        %CheckArray(2*x-3,y) = CheckArray(2*x-3,y) + 1;
        Fgrid(2*x-2,y,:,:) = Matrix_FourPointAvarage(squeeze(grid(x-1,n-1,:,:)) , squeeze(grid(x,n-1,:,:)) , squeeze(grid(x+1,n-1,:,:)) , squeeze(grid(x+2,n-1,:,:)),omega); 
        %CheckArray(2*x-2,y) = CheckArray(2*x-2,y) + 1;
end
Lastline(1,y,:,:) = squeeze(grid(m,n,:,:));
% Here - 'x' fixed
x = 2*m-5;
for y=2:(n-2)
        Fgrid(x,2*y-3,:,:) = squeeze(grid(m-1,y,:,:));
        %CheckArray(x,2*y-3) = CheckArray(x,2*y-3) + 1;
        Fgrid(x,2*y-2,:,:) = Matrix_FourPointAvarage(squeeze(grid(m-1,y-1,:,:)) , squeeze(grid(m-1,y,:,:)) , squeeze(grid(m-1,y+1,:,:)) , squeeze(grid(m-1,y+2,:,:)),omega); 
        %CheckArray(x,2*y-2) = CheckArray(x,2*y-2) + 1;
end
%-----------------------------------------------------------

%-----------------------------------------------------------
% Second main loops : vertices of face center 
for x=3:(m-3)
    for y=3:(n-3)
        % X direction
        m1 = Matrix_FourPointAvarage(squeeze(Fgrid(2*x-5,2*y-2,:,:)) , squeeze(Fgrid(2*x-3,2*y-2,:,:)) , squeeze(Fgrid(2*x-1,2*y-2,:,:)) , squeeze(Fgrid(2*x+1,2*y-2,:,:)),omega);
        % Y direction
        m2 = Matrix_FourPointAvarage(squeeze(Fgrid(2*x-2,2*y-5,:,:)) , squeeze(Fgrid(2*x-2,2*y-3,:,:)) , squeeze(Fgrid(2*x-2,2*y-1,:,:)) , squeeze(Fgrid(2*x-2,2*y+1,:,:)),omega);
        % in scalars or commutative space m1=m1. Here is not!
        Fgrid(2*x-2,2*y-2,:,:) = GM2(m1,m2,.5);
        %CheckArray(2*x-2,2*y-2) = CheckArray(2*x-2,2*y-2) + 1;
    end
end
%-----------------------------------------------------------
% Bounderies handeling - vertices of face center
%-----------------------------------------------------------
% Line 2:
x=2;
for y=2:(n-3)
    % X direction
    m1 = Matrix_FourPointAvarage(squeeze(Firstline(1,2*y-2,:,:)) , squeeze(Fgrid(2*x-3,2*y-2,:,:)) , squeeze(Fgrid(2*x-1,2*y-2,:,:)) , squeeze(Fgrid(2*x+1,2*y-2,:,:)),omega);
    % Y direction
    m2 = Matrix_FourPointAvarage(squeeze(Firstrow(1,2*y-2,:,:)) , squeeze(Fgrid(2*x-2,2*y-3,:,:)) , squeeze(Fgrid(2*x-2,2*y-1,:,:)) , squeeze(Fgrid(2*x-2,2*y+1,:,:)),omega);
    % in scalars or commutative space m1=m1. Here is not!
    Fgrid(2*x-2,2*y-2,:,:) = GM2(m1,m2,.5);
    %CheckArray(2*x-2,2*y-2) = CheckArray(2*x-2,2*y-2) + 1;
end

% HERE iS last place using last row
%-----------------------------------------------------------
% Line 2 before end:
x=m-2;
for y=2:(n-3)
    % X direction
    m1 = Matrix_FourPointAvarage(squeeze(Fgrid(2*x-5,2*y-2,:,:)) , squeeze(Fgrid(2*x-3,2*y-2,:,:)) , squeeze(Fgrid(2*x-1,2*y-2,:,:)) , squeeze(Lastline(1,2*y-2,:,:)),omega);
    % Y direction
    m2 = Matrix_FourPointAvarage(squeeze(Firstrow(1,2*y-2,:,:)) , squeeze(Fgrid(2*x-2,2*y-3,:,:)) , squeeze(Fgrid(2*x-2,2*y-1,:,:)) , squeeze(Fgrid(2*x-2,2*y+1,:,:)),omega);
    % in scalars or commutative space m1=m1. Here is not!
    Fgrid(2*x-2,2*y-2,:,:) = GM2(m1,m2,.5);
    %CheckArray(2*x-2,2*y-2) = CheckArray(2*x-2,2*y-2) + 1;
end

% HERE iS last place using last row
%-----------------------------------------------------------
% Row 2:
y=2;
for x=3:(m-3)
    % X direction
    m1 = Matrix_FourPointAvarage(squeeze(Firstline(1,2*y-2,:,:)) , squeeze(Fgrid(2*x-3,2*y-2,:,:)) , squeeze(Fgrid(2*x-1,2*y-2,:,:)) , squeeze(Fgrid(2*x+1,2*y-2,:,:)),omega);
    % Y direction
    m2 = Matrix_FourPointAvarage(squeeze(Firstrow(1,2*y-2,:,:)) , squeeze(Fgrid(2*x-2,2*y-3,:,:)) , squeeze(Fgrid(2*x-2,2*y-1,:,:)) , squeeze(Fgrid(2*x-2,2*y+1,:,:)),omega);
    % in scalars or commutative space m1=m1. Here is not!
    Fgrid(2*x-2,2*y-2,:,:) = GM2(m1,m2,.5);
    %CheckArray(2*x-2,2*y-2) = CheckArray(2*x-2,2*y-2) + 1;
end
%-----------------------------------------------------------
% Row 2 before end:
y=n-2;
for x=3:(n-3)
    % X direction
    m1 = Matrix_FourPointAvarage(squeeze(Fgrid(2*x-5,2*y-2,:,:)) , squeeze(Fgrid(2*x-3,2*y-2,:,:)) , squeeze(Fgrid(2*x-1,2*y-2,:,:)) , squeeze(Lastline(1,2*y-2,:,:)),omega);
    % Y direction
    m2 = Matrix_FourPointAvarage(squeeze(Fgrid(2*x-2,2*y-5,:,:)) , squeeze(Fgrid(2*x-2,2*y-3,:,:)) , squeeze(Fgrid(2*x-2,2*y-1,:,:)) , squeeze(Lastrow(1,2*x-2,:,:)),omega);
    % in scalars or commutative space m1=m1. Here is not!
    Fgrid(2*x-2,2*y-2,:,:) = GM2(m1,m2,.5);
    %CheckArray(2*x-2,2*y-2) = CheckArray(2*x-2,2*y-2) + 1;
end
%-----------------------------------------------------------
%Last two cases:
%-----------------------------------------------------------
x = 2;
y = n-2;
    % X direction
    m1 = Matrix_FourPointAvarage(squeeze(Firstline(1,2*y-2,:,:)) , squeeze(Fgrid(2*x-3,2*y-2,:,:)) , squeeze(Fgrid(2*x-1,2*y-2,:,:)) , squeeze(Fgrid(2*x+1,2*y-2,:,:)),omega);
    % Y direction
    m2 = Matrix_FourPointAvarage(squeeze(Fgrid(2*x-2,2*y-5,:,:)) , squeeze(Fgrid(2*x-2,2*y-3,:,:)) , squeeze(Fgrid(2*x-2,2*y-1,:,:)) , squeeze(Lastrow(1,2*x-2,:,:)),omega);
    % in scalars or commutative space m1=m1. Here is not!
    Fgrid(2*x-2,2*y-2,:,:) = GM2(m1,m2,.5);
    %CheckArray(2*x-2,2*y-2) = CheckArray(2*x-2,2*y-2) + 1;
%-----------------------------------------------------------
x = m-2;
y = n-2;
    % X direction
    m1 = Matrix_FourPointAvarage(squeeze(Fgrid(2*x-5,2*y-2,:,:)) , squeeze(Fgrid(2*x-3,2*y-2,:,:)) , squeeze(Fgrid(2*x-1,2*y-2,:,:)) , squeeze(Lastline(1,2*y-2,:,:)),omega);
    % Y direction
    m2 = Matrix_FourPointAvarage(squeeze(Lastline(1,2*y-5,:,:)) , squeeze(Fgrid(2*x-2,2*y-3,:,:)) , squeeze(Fgrid(2*x-2,2*y-1,:,:)) , squeeze(Lastrow(1,2*x-2,:,:)),omega);
    % in scalars or commutative space m1=m1. Here is not!
    Fgrid(2*x-2,2*y-2,:,:) = GM2(m1,m2,.5);
    %CheckArray(2*x-2,2*y-2) = CheckArray(2*x-2,2*y-2) + 1;
%-----------------------------------------------------------



end

